Univariate and Bivariate Analyses

source("~/Documents/GitHub/VivID_Epi/analyses/00-functions.R") 
library(tidyverse)
library(sf)
library(srvyr) #wrap the survey package in dplyr syntax
#......................
# Import Data
#......................
load("~/Documents/GitHub/VivID_Epi/data/vividepi_raw.rda")
dt <- merge_pr_plsmdm_gemtdt(pr = arpr, plsmdm = panplasmpcrres, ge = ge)

Pf versus Pv Prevalence

Steve has asked me to produce maps and basic correlation plots of Pf versus Pv

prev_summarizer <- function(data, maplvl, plsmdmspec, sfobj){
  
  
  # catch error
  # if( !( any(colnames(x) %in% c("hv001")) & any(colnames(x) %in% c("hv005")) ) ){
  #   stop("Must include cluster and cluster weight. need to check if you are using pr or hiv sample weights")
  # }
  # rlang
  maplvl <- enquo(maplvl)
  plsmdmspec <- enquo(plsmdmspec)
  
  # clusters are weighted (each individual has same weight in cluster)
  ret <- data %>% 
    srvyr::as_survey_design(ids = hv001, weights = hv005) %>% 
    dplyr::group_by(!!maplvl) %>% 
    dplyr::summarise(plsmdn = srvyr::survey_total(hv001), 
                     plsmd = srvyr::survey_mean(!!plsmdmspec, na.rm = T, vartype = c("se", "ci"), level = 0.95)) %>% 
    dplyr::left_join(., sfobj) # attach spatial data, let R figure out the common var
  # return
  return(ret)
}



pfldhprov <- prev_summarizer(data = dt, maplvl = adm1name, plsmdmspec = pfldh, sfobj = DRCprov) %>% 
  dplyr::mutate(plsmdmspec = "pfldh", maplvl = "adm1name")
pv18sprov <- prev_summarizer(data = dt, maplvl = adm1name, plsmdmspec = pv18s, sfobj = DRCprov)  %>% 
  dplyr::mutate(plsmdmspec = "pv18s", maplvl = "adm1name")

pfldhclust <- prev_summarizer(data = dt, maplvl = hv001, plsmdmspec = pfldh, sfobj = ge) %>% 
  dplyr::mutate(plsmdmspec = "pfldh", maplvl = "hv001")
pv18sclust <- prev_summarizer(data = dt, maplvl = hv001, plsmdmspec = pv18s, sfobj = ge) %>% 
  dplyr::mutate(plsmdmspec = "pv18s", maplvl = "hv001")


mp <- dplyr::bind_rows(pfldhprov, pv18sprov, pfldhclust, pv18sclust) %>% 
  dplyr::group_by(plsmdmspec, maplvl) %>% 
  tidyr::nest()

mp$data <- sapply(list(pfldhprov, pv18sprov, pfldhclust, pv18sclust), function(x) return(x))

mapplotter <- function(data, maplvl, plsmdmspec){
  
  
  ret <- ggplot() + 
    geom_sf(data = DRCprov) +
    ggtitle(paste(plsmdmspec)) +
    theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 13), 
          legend.position = "bottom",
          legend.title = element_text(famil = "Arial", face = "bold", vjust = 0.85, size = 12),
          legend.text = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
          axis.text = element_blank(),
          axis.line = element_blank(),
          panel.background = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent"),
          panel.grid = element_blank(),
          panel.border = element_blank())
  
  if(maplvl == "adm1name"){
    
    ret <- ret + geom_sf(data = data, aes(fill = plsmd)) +
      scale_fill_gradient2("Prevalence", low = "#0000FF", mid = "#FFEC00", high = "#FF0000") + 
      coord_sf(datum=NA)  # to get rid of gridlines
    
  } else if(maplvl == "hv001"){
    
    ret <- ret + geom_sf(data = data, aes(fill = plsmd, colour = plsmd, size = plsmdn), alpha = 0.8) +
      scale_color_gradient2("Prevalence", low = "#0000FF", mid = "#FFEC00", high = "#FF0000") + 
      scale_size(guide = 'none') +  scale_fill_continuous(guide = 'none') +
      coord_sf(datum=NA)  # to get rid of gridlines
    
  } else {
    stop("maplvl is not in the options for this function")
  }
  
  
  return(ret)
  
}

prevmaps <- pmap(mp, mapplotter)

Maps

Pfalciparum Prevalence

gridExtra::grid.arrange(prevmaps[[1]], 
                        prevmaps[[3]],
                        ncol=2, top=grid::textGrob("Pf Prevalence by Province and Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) 

Pvivax Prevalence

gridExtra::grid.arrange(prevmaps[[2]], 
                        prevmaps[[4]],
                        ncol=2, top=grid::textGrob("Pf Prevalence by Province and Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) 

Pvivax and Pfalciparum Prevalence

gridExtra::grid.arrange(prevmaps[[1]], 
                        prevmaps[[2]],
                        prevmaps[[3]],
                        prevmaps[[4]],
                        ncol=2, nrow=2, top=grid::textGrob("Prevalence by Species and Province/Cluster in CD2013 DHS", gp=grid::gpar(fontsize=15, fontfamily = "Arial", fontface = "bold"))) 

Pv versus Pf Correlations

By Province

Points are scaled by size of province (and color of province). Loess curve in red.

pfprov <- mp$data[[1]]
colnames(pfprov)[4] <- "Pf_prev"

pvprov <- mp$data[[2]]
colnames(pvprov)[4] <- "Pv_prev"
  

plotObj <- dplyr::bind_cols(pfprov, pvprov) %>% 
  ggplot(aes(x=Pf_prev, y=Pv_prev, size = plsmdn)) +
  geom_point(aes(colour = adm1name)) +
  geom_smooth(method="loess", se=F, colour = "red") +
      ggtitle("Pfalciparum versus Pvivax Province Prevalence") + ylab("Pv Prevalence") + xlab("Pf Prevalence") +
    theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 18), 
          legend.position = "none",
          axis.text.x = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
          axis.text.y = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, size = 10),
          axis.ticks = element_blank(),
          axis.line = element_line(colour = "black"),
          panel.background = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent"),
          panel.grid = element_blank(),
          panel.border = element_blank())

plotly::ggplotly(plotObj)

By Cluster

Points are scaled by size of province (and color of province). Loess curve in red. Cluster filtered for both 0s.

pfprov <- mp$data[[3]]
colnames(pfprov)[4] <- "Pf_prev"

pvprov <- mp$data[[4]]
colnames(pvprov)[4] <- "Pv_prev"
  

plotObj <- dplyr::bind_cols(pfprov, pvprov) %>% 
  dplyr::filter(!(Pf_prev == 0 & Pv_prev == 0)) %>% 
  ggplot(aes(x=Pf_prev, y=Pv_prev, size = plsmdn)) +
  geom_point(aes(colour = hv001)) +
  geom_smooth(method="loess", se=F, colour = "red") +
      ggtitle("Pfalciparum versus Pvivax Province Prevalence") + ylab("Pv Prevalence") + xlab("Pf Prevalence") +
    theme(plot.title = element_text(famil = "Arial", face = "bold", hjust = 0.5, size = 18), 
          legend.position = "none",
          axis.text.x = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, angle = 90, size = 10),
          axis.text.y = element_text(famil = "Arial", hjust = 0.5, vjust = 0.5, size = 10),
          axis.ticks = element_blank(),
          axis.line = element_line(colour = "black"),
          panel.background = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent"),
          panel.grid = element_blank(),
          panel.border = element_blank())

plotly::ggplotly(plotObj)